# Current, as of Thursday, May 28, 2020.
# rm(list = ls())
# setwd("~/Dropbox/Base Dropbox/List Experiment Paper (PS)/Conditional Accept/Replication Files (Carmines-Schmidt PS)/2012 CCAP/Models (full sample)")
library(foreign)
library(car)
library(readstata13)
library(zeligverse)
library(ggplot2)
#-----
# Input the data.
our.data = read.dta13("CarminesSchmidtPS-replication-2012.dta", generate.factors = TRUE)
names(our.data)
#-----
# Subsets, based on assignment to the "baseline" condition or each sensitive-item condition. 

# "A candidate who is AFRICAN AMERICAN" (treatment) or baseline condition (dropping black respondents from the sample). 
afroam.experiment = subset(our.data, treatment=="African Americans" | treatment=="Control")
afroam.experiment = subset(afroam.experiment, black!=1)

# "A candidate who is MUSLIM" (treatment) or baseline condition. 
muslims.experiment = subset(our.data, treatment=="Muslims" | treatment=="Control")

# "A candidate who is GAY OR HOMOSEXUAL" (treatment) or baseline condition. 
gays.experiment = subset(our.data, treatment=="Gays" | treatment=="Control")

# "A candidate who is MORMON" (treatment) or baseline condition. 
mormons.experiment = subset(our.data, treatment=="Mormons" | treatment=="Control")
#-----
# Treatment effect: AFRICAN AMERICAN list experiment. 
set.seed(47408) 
our.model = lm(item_count ~ afroam_binary, data = afroam.experiment)
summary(our.model)
# Use Zelig to model the difference in means based on treatment assignment.. 
z.out = zelig(item_count ~ afroam_binary,  
              model = "normal", data = afroam.experiment)
z.out
# Not assigned to the sensitive-item condition. 
x.low = setx(z.out, 
             afroam_binary = 0)
# Assigned to the sensitive-item condition. 
x.high = setx(z.out, 
              afroam_binary = 1)
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and 90-percent confidence intervals representing the treatment effect. 
LB.afroam = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.afroam = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.afroam = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))
#-----
# Treatment effect: GAY OR HOMOSEXUAL list experiment. 
set.seed(47408)
our.model = lm(item_count ~ gays_binary, data = gays.experiment)
summary(our.model)
# Model item count as a function of assignment to the sensitive-item condition.
z.out = zelig(item_count ~ gays_binary,  
              model = "normal", data = gays.experiment)
z.out
# Not assigned to the sensitive-item condition. 
x.low = setx(z.out, 
             gays_binary = 0)
# Assigned to the sensitive-item condition. 
x.high = setx(z.out, 
              gays_binary = 1)
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and 90-percent confidence intervals representing the treatment effect. 
LB.gays = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.gays = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.gays = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))
#-----
# Treatment effect: MUSLIM list experiment. 
set.seed(47408)
our.model = lm(item_count ~ muslim_binary, data = muslims.experiment)
summary(our.model)
# Model item count as a function of assignment to the sensitive-item condition.
z.out = zelig(item_count ~ muslim_binary,  
              model = "normal", data = muslims.experiment)
z.out
# Not assigned to the sensitive-item condition. 
x.low = setx(z.out, 
             muslim_binary = 0)
# Assigned to the sensitive-item condition. 
x.high = setx(z.out, 
              muslim_binary = 1)
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and 90-percent confidence intervals representing the treatment effect. 
LB.muslims = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.muslims = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.muslims = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))
#-----
# Treatment effect: MORMON list experiment. 
set.seed(47408)
our.model = lm(item_count ~ mormon_binary, data = mormons.experiment)
summary(our.model)
# Model item count as a function of assignment to the sensitive-item condition.
z.out = zelig(item_count ~ mormon_binary,  
              model = "normal", data = mormons.experiment)
z.out
# Not assigned to the sensitive-item condition. 
x.low = setx(z.out, 
             mormon_binary = 0)
# Assigned to the sensitive-item condition. 
x.high = setx(z.out, 
              mormon_binary = 1)
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals representing the treatment effect. 
LB.mormon = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.mormon = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.mormon = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))
#-----
# Point estimates.
# round(cbind(PE.afroam, PE.gays, PE.muslims, PE.mormon), digits = 2)
# Lower bounds.
# round(cbind(LB.afroam, LB.gays, LB.muslims, LB.mormon), digits = 2)
# Upper bounds.
# round(cbind(UB.afroam, UB.gays, UB.muslims, UB.mormon), digits = 2)
#-----
# Treatment effect: AFRICAN AMERICAN list experiment, including control variables. 
set.seed(47408)
our.model = lm(item_count ~ afroam_binary + pid_7 + ideology + female + white + 
                 education + church + rr + obama_favor + interest + income, 
               data = afroam.experiment)
summary(our.model)
# Model item count as a function of assignment to the sensitive-item condition.
z.out = zelig(item_count ~ afroam_binary + pid_7 + ideology + female + white + 
                education + church + rr + obama_favor + interest + income, 
              model = "normal", data = afroam.experiment)
z.out
# Not assigned to the sensitive-item condition. 
x.low = setx(z.out, 
             afroam_binary = 0,
             pid_7 = mean(afroam.experiment$pid_7, na.rm=TRUE),
             ideology = mean(afroam.experiment$ideology, na.rm=TRUE),
             female = median(afroam.experiment$female, na.rm=TRUE),
             white = median(afroam.experiment$white, na.rm=TRUE),
             education = median(afroam.experiment$education, na.rm=TRUE),
             church = median(afroam.experiment$church, na.rm=TRUE),
             rr = mean(afroam.experiment$rr, na.rm=TRUE),
             obama_favor = median(afroam.experiment$obama_favor, na.rm=TRUE),
             interest = median(afroam.experiment$interest, na.rm=TRUE),
             income = median(afroam.experiment$income, na.rm=TRUE))
# Assigned to the sensitive-item condition. 
x.high = setx(z.out, 
              afroam_binary = 1,
              pid_7 = mean(afroam.experiment$pid_7, na.rm=TRUE),
              ideology = mean(afroam.experiment$ideology, na.rm=TRUE),
              female = median(afroam.experiment$female, na.rm=TRUE),
              white = median(afroam.experiment$white, na.rm=TRUE),
              education = median(afroam.experiment$education, na.rm=TRUE),
              church = median(afroam.experiment$church, na.rm=TRUE),
              rr = mean(afroam.experiment$rr, na.rm=TRUE),
              obama_favor = median(afroam.experiment$obama_favor, na.rm=TRUE),
              interest = median(afroam.experiment$interest, na.rm=TRUE),
              income = median(afroam.experiment$income, na.rm=TRUE))
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals representing the treatment effect. 
LB.afroam = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.afroam = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.afroam = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))
#-----
# Treatment effect: GAY OR HOMOSEXUAL list experiment, including control variables. 
set.seed(47408)
our.model = lm(item_count ~ gays_binary + pid_7 + ideology + female + white + 
                 black + education + church + rr + obama_favor + interest + income, 
               data = gays.experiment)
summary(our.model)
# Model item count as a function of assignment to the sensitive-item condition.
z.out = zelig(item_count ~ gays_binary + pid_7 + ideology + female + white + 
                black + education + church + rr + obama_favor + interest + income,  
              model = "normal", data = gays.experiment)
z.out
# Not assigned to the sensitive-item condition. 
x.low = setx(z.out, 
             gays_binary = 0,
             pid_7 = mean(gays.experiment$pid_7, na.rm=TRUE),
             ideology = mean(gays.experiment$ideology, na.rm=TRUE),
             female = median(gays.experiment$female, na.rm=TRUE),
             white = median(gays.experiment$white, na.rm=TRUE),
             black = median(gays.experiment$black, na.rm=TRUE),
             education = median(gays.experiment$education, na.rm=TRUE),
             church = median(gays.experiment$church, na.rm=TRUE),
             rr = mean(gays.experiment$rr, na.rm=TRUE),
             obama_favor = median(gays.experiment$obama_favor, na.rm=TRUE),
             interest = median(gays.experiment$interest, na.rm=TRUE),
             income = median(gays.experiment$income, na.rm=TRUE))
# Assigned to the sensitive-item condition. 
x.high = setx(z.out, 
              gays_binary = 1,
              pid_7 = mean(gays.experiment$pid_7, na.rm=TRUE),
              ideology = mean(gays.experiment$ideology, na.rm=TRUE),
              female = median(gays.experiment$female, na.rm=TRUE),
              white = median(gays.experiment$white, na.rm=TRUE),
              black = median(gays.experiment$black, na.rm=TRUE),
              education = median(gays.experiment$education, na.rm=TRUE),
              church = median(gays.experiment$church, na.rm=TRUE),
              rr = mean(gays.experiment$rr, na.rm=TRUE),
              obama_favor = median(gays.experiment$obama_favor, na.rm=TRUE),
              interest = median(gays.experiment$interest, na.rm=TRUE),
              income = median(gays.experiment$income, na.rm=TRUE))
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals representing the treatment effect. 
LB.gays = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.gays = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.gays = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))
#-----
# Treatment effect: MUSLIM list experiment, including control variables. 
set.seed(47408)
our.model = lm(item_count ~ muslim_binary + pid_7 + ideology + female + white + 
                 black + education + church + rr + obama_favor + interest + income, 
               data = muslims.experiment)
summary(our.model)
# Model item count as a function of assignment to the sensitive-item condition.
z.out = zelig(item_count ~ muslim_binary + pid_7 + ideology + female + white +
                black + education + church + rr + obama_favor + interest + income,  
              model = "normal", data = muslims.experiment)
z.out
# Not assigned to the sensitive-item condition. 
x.low = setx(z.out, 
             muslim_binary = 0,
             pid_7 = mean(muslims.experiment$pid_7, na.rm=TRUE),
             ideology = mean(muslims.experiment$ideology, na.rm=TRUE),
             female = median(muslims.experiment$female, na.rm=TRUE),
             white = median(muslims.experiment$white, na.rm=TRUE),
             black = median(muslims.experiment$black, na.rm=TRUE),
             education = median(muslims.experiment$education, na.rm=TRUE),
             church = median(muslims.experiment$church, na.rm=TRUE),
             rr = mean(muslims.experiment$rr, na.rm=TRUE),
             obama_favor = median(muslims.experiment$obama_favor, na.rm=TRUE),
             interest = median(muslims.experiment$interest, na.rm=TRUE),
             income = median(muslims.experiment$income, na.rm=TRUE))
# Assigned to the sensitive-item condition. 
x.high = setx(z.out, 
              muslim_binary = 1,
              pid_7 = mean(muslims.experiment$pid_7, na.rm=TRUE),
              ideology = mean(muslims.experiment$ideology, na.rm=TRUE),
              female = median(muslims.experiment$female, na.rm=TRUE),
              white = median(muslims.experiment$white, na.rm=TRUE),
              black = median(muslims.experiment$black, na.rm=TRUE),
              education = median(muslims.experiment$education, na.rm=TRUE),
              church = median(muslims.experiment$church, na.rm=TRUE),
              rr = mean(muslims.experiment$rr, na.rm=TRUE),
              obama_favor = median(muslims.experiment$obama_favor, na.rm=TRUE),
              interest = median(muslims.experiment$interest, na.rm=TRUE),
              income = median(muslims.experiment$income, na.rm=TRUE))
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals representing the treatment effect. 
LB.muslims = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.muslims = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.muslims = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))
#-----
# Treatment effect: MORMON list experiment, including control variables. 
set.seed(47408)
our.model = lm(item_count ~ mormon_binary + pid_7 + ideology + female + white + 
                 black + education + church + rr + obama_favor + interest + income + 
                 romney_favor, 
               data = mormons.experiment)
summary(our.model)
# Model item count as a function of assignment to the sensitive-item condition.
z.out = zelig(item_count ~ mormon_binary + pid_7 + ideology + female + white + 
                black + education + church + rr + obama_favor + interest + income + 
                romney_favor, 
              model = "normal", data = mormons.experiment)
z.out
# Not assigned to the sensitive-item condition. 
x.low = setx(z.out, 
             mormon_binary = 0,
             pid_7 = mean(mormons.experiment$pid_7, na.rm=TRUE),
             ideology = mean(mormons.experiment$ideology, na.rm=TRUE),
             female = median(mormons.experiment$female, na.rm=TRUE),
             white = median(mormons.experiment$white, na.rm=TRUE),
             black = median(mormons.experiment$black, na.rm=TRUE),
             education = median(mormons.experiment$education, na.rm=TRUE),
             church = median(mormons.experiment$church, na.rm=TRUE),
             rr = mean(mormons.experiment$rr, na.rm=TRUE),
             obama_favor = median(mormons.experiment$obama_favor, na.rm=TRUE),
             interest = median(mormons.experiment$interest, na.rm=TRUE),
             income = median(mormons.experiment$income, na.rm=TRUE),
             romney_favor = median(mormons.experiment$romney_favor, na.rm=TRUE))
# Assigned to the sensitive-item condition. 
x.high = setx(z.out, 
              mormon_binary = 1,
              pid_7 = mean(mormons.experiment$pid_7, na.rm=TRUE),
              ideology = mean(mormons.experiment$ideology, na.rm=TRUE),
              female = median(mormons.experiment$female, na.rm=TRUE),
              white = median(mormons.experiment$white, na.rm=TRUE),
              black = median(mormons.experiment$black, na.rm=TRUE),
              education = median(mormons.experiment$education, na.rm=TRUE),
              church = median(mormons.experiment$church, na.rm=TRUE),
              rr = mean(mormons.experiment$rr, na.rm=TRUE),
              obama_favor = median(mormons.experiment$obama_favor, na.rm=TRUE),
              interest = median(mormons.experiment$interest, na.rm=TRUE),
              income = median(mormons.experiment$income, na.rm=TRUE),
              romney_favor = median(mormons.experiment$romney_favor, na.rm=TRUE))
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals representing the treatment effect. 
LB.mormon = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.mormon = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.mormon = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))
#-----
# Point estimates.
# round(cbind(PE.afroam, PE.gays, PE.muslims, PE.mormon), digits = 2)
# Lower bounds.
# round(cbind(LB.afroam, LB.gays, LB.muslims, LB.mormon), digits = 2)
# Upper bounds.
# round(cbind(UB.afroam, UB.gays, UB.muslims, UB.mormon), digits = 2)


